function [f,g]=LLHs(beta,X,TC_idx,q0,q1,sqw)
KX=size(X,2); L=length(q0);
save tmpdos.mat beta
% Build vectors of expressions containing parameters
delta=X*beta;  ed=exp(delta); % exponential of delta
if sum(isinf(ed)) > 0 && sum(beta(1:3)>0) >1
   fprintf('pause since delta is > 900 and 2 of 3 beta are positive');
    pause
end
sed=my_accumarray(TC_idx,ed); %sum of ed(j) within group TC
pr=ed./(1+sed); 
lpr=log(pr); q1lpr=q1.*lpr; q1lpr(q0)=0;
% Calculate the LLH
f1=sqw'*q1lpr;
f= - f1; %minimize -LLH
if sum(isnan(f)) > 0
   fprintf('pause since f is nan');
    pause
end
if sum(isinf(f)) > 0
   fprintf('pause since f is inf');
    pause
end
if sum(~isreal(f)) > 0
   fprintf('pause since f is complex');
    pause
end
if nargout==1, return; end
% Calculate gradient parts
prX=pr.*X; com0=my_accumarray(TC_idx,prX);
g=sqw'*(X-com0);
g= - g; %minimize -LLH
if (sum(sum(isnan(g)))+sum(sum(isinf(g)))+sum(sum(~isreal(g)))) > 0
    find(isnan(g))
    find(isinf(g))
    find(~isreal(g))
   fprintf('pause since g is nan or inf or complex');
    pause
end
end
